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Abstract 

This  paper  discusses  the  application  of  the  three 
optimization  methods  —  steepest  descent,  second  varia¬ 
tion,  and  generalized  Newton-Raphson  —  to  the  problem  of 
minimum  time,  low  thrust,  circle- to-circle  transfer.  De¬ 
tails  of  computational  techniques  that  have  proved  suc¬ 
cessful  in  practice  are  presented.  The  number  of  itera¬ 
tion  cycles  and  the  time  used  by  the  computer  are  given 
for  each  method. 
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Introduction  and  Statement  of  Problem 


This  paper  describes  the  application  of  the  three 
basic  optimization  schemes  covered  in  Part  I:  Discussion 
to  a  specific  problem.  This  problem  is  that  of  minimiz¬ 
ing  the  transfer  time  of  a  low-thrust  ion  rocket  between 
the  orbits  of  Earth  and  Mars.  Applications  of  the  gradi¬ 
ent  method  and  of  the  second  variation  method  to  this 
problem  have  been  reported  in  QJ  ,  (2),  and  (4),  and 
therefore  only  a  brief  sketch  of  interesting  variations 
in  these  computational  procedures  and  of  the  results  will 
be  included  here. 

To  simplify  the  problem  as  much  as  possible  the 
rocket's  thrust  level  was  assumed  constant,  and  thus  the 
single  control  variable  is  the  thrust  direction.  Further, 
the  orbits  of  Earth  and  Mars  were  assumed  to  be  circular 
and  coplanar,  and  the  gravitational  attractions  of  the 
two  planets  on  the  vehicle  were  neglected.  The  following 
system  parameters  for  the  low-thrust  vehicle  were  adopted 
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from  (1)  . 


Initial  Mass,  ihq 
Specific  Impulse 
Propellant  Consumption,  m, 
Thrust,  T, 

Tlirust/ Initial  Weight 


46.58  slugs 
5700  sec 

-6.937  x  10  7  slugs/sec 
0.127  lb 
0.9  x  10“4 


The  equations  of  motion  were  given  by: 
Radial  Velocity 


=  £w  - 


=  u 


Radial  Acceleration 

f(2)  vl 

u  =  f  J  -  —  ■ 
r 


Ji_  ,  I-.sin__ 
r2  iDq  +  mt 


Circumferential  Acceleration 


(3)  =  - 


it\  +  T  cos 

r  niQ  +  mt 


(1) 

(2) 

(3) 


where  u  and  v  are  the  radial  and  circumferential 
velocities  respectively;  r  is  the  radius;  and  is 

the  thrust  direction  angle  measured  from  the  local  hori¬ 
zontal.  All  the  initial  and  final  values  of  the  state 
variables  were  specified,  and  the  quantity  to  be  mini¬ 
mized  was  t£,  the  final  time. 


The  Three  Methods  as  Applied  to  the 
Sample  Problem 


Gradient  Method 


As  stated  in  Part  I,  the  gradient  methods  utilize  a 
penalty  function  which  for  this  problem  becomes 


(4) 
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As  was  also  stated  in  Part  I,  the  criterion  for  the  ter 
mination  of  a  trajectory  is 


P  = 


V 

L 

i=l 


1  +  )  k.  (x . 

L  11, 


x. )  f 
i  f 


=  0 


Utilizing 


A 


f  l 


=  k.  (x.  -  x.) 

ox.  l  ir  l ' 


we  have 


1  +  l  A t  f =  0 

i=l  r 


(5) 


(6) 


(7) 


This  equation  scales  the  multipliers. 

Since  by  the  "multiplier  rule"  all  the  final  A's 
cannot  vanish,  the  extremal,  that  minimizes  P1  must  de¬ 
viate  from  the  desired  values  x^ .  We  will  present  two 
methods  of  reducing  the  deviation  to  within  tolerable 
limits . 

We  will  assume  that  the  errors  are  small  enough  so 
that  the  extremal  that  minimizes  P1  is  close  to  the 
extremal  that  passes  through  the  desired  end-point.  By 
close  we  mean  that  the  two  extremals  have  approximately 
the  same  final  A's.  If  we  increase  the  k' s  and  re¬ 
converge,  this  new  extremal  will  also  have  approximately 
the  same  final  A's.  It  then  follows  from  Eq.  (6)  that 
the  deviations  will  be  reduced  in  proportion  to  the  in¬ 
crease  in  the  k's. 

Thus  the  first  method  of  reducing  the  errors  is  to 
"build  up  the  k's."  This  method  was  often  unsuccessful 
however,  due  to  inexact  interpolation  for  the  final 
values  and  other  numerical  errors  that  prevented  conver¬ 
gence.  Therefore,  the  method  of  "shifting  boundary 
values"  was  devised  and  has  proved  invariably  successful. 

After  convergence,  the  boundary  values  are  rede¬ 
fined  by  subtracting  Ax^  =  xj.  -  x^  from  x^.  We  then 

apply  the  gradient  process  to  obtain  the  minimum  of 
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3 

P  =  t  +  y  k^(x^  -  +  Ax^)^  (8) 

i=l  f 

If  the  final  A' s  were  unchanged,  the  extremal  would 
pass  directly  through  the  desired  end-point.  In  prac¬ 
tice,  several  boundary  shifts  may  be  necessary. 

Applying  a  gradient  scheme  incorporating  the  bound¬ 
ary  shift  technique  to  the  above  problem  of  minimizing 
the  time  of  low-thrust  Earth  to  Mars  trajectories,  it  was 
found  that  attainment  of  the  minimum  t  of  193  days 

Vvith  convergence  to  within  0.17o  of  the  boundary  values 
was  possible  in  5  shifts  with  an  average  of  20  descents 
per  shift.  As  programmed,  a  descent  was  a  one  dimen¬ 
sional  search  for  a  minimum  of  P1  versus  control  varia¬ 
tions  calculated  from  the  adjoint  equations,  and  it  re¬ 
quired  the  computation  of  one  adjoint  solution  and  at 
most  four  trajectories.  Figure  1  shows  the  initial, 
some  intermediate,  and  the  final  Q  time-histories  pre¬ 
ceding  the  first  boundary  shift.  As  indicated  in  the 
figure,  the  first  boundary  shift  took  place  after  25  de¬ 
scents.  The  converged  curve,  labeled  F,  is  also  in¬ 
cluded. 

Second  Variation  Method 

As  explained  in  Part  I,  the  first  stage  of  the  sec¬ 
ond  variation  method  also  depends  on  a  penalty  function, 
P' .  However,  in  this  case  P1  is  expanded  to  second 
order  terms,  and  among  these  are  terms  in  6tf,  the 
change  in  the  final  time.  Because  these  terms  greatly 
complicate  an  already  lengthy  numerical  calculation, 

6tf  was  set  equal  to  zero  in  the  expansion  of  P1  when 

the  second  variation  method  was  coded  for  the  computer. 
This  implies  that  the  control  variable  increment  60  (t) 
was  then  calculated  to  obtain  the  maximum  reduction  in 
P  at  the  nominal  final  time  tf  instead  of  tf  +  5tf. 

Of  course  when  the  trajectory  with  control  variable 
v  +  60  was  computed,  it  was  not  terminated  at  time  tf 
but  rather  at  that  point  of  the  trajectory  with  minimum 
P1,  as  is  done  in  the  gradient  programs.  With  respect 
to  over-all  computational  time  this  technique  represents 
a  compromise  between  a  true  second  variation  calculation 
and  additional  programming  complexity.  For  this  particu¬ 
lar  problem  it  was  advantageous  to  treat  the  problem  as  a 
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fixed  time  problem  when  computing  68  using  penalty 
functions . 

The  initial  0 (t)  function  in  this  case  corre¬ 
sponds  to  constant  circumferential  thrust.  This  resulted 
in  terminal  boundary  value  errors  that  averaged  20  per 
cent.  After  6  descent  cycles,  using  the  penalty  function 
procedure,  the  terminal  errors  averaged  3  per  cent  with 
the  transfer  time  at.  180  days.  After  5  additional  cycles 
of  the  refinement  process,  the  average  boundary  value 
error  was  reduced  to  0.05  per  cent  and  the  transfer  time 
had  reached  its  minimum  of  193  days.  The  over-all 
IBM  7094  computer  time  was  2  minutes,  representing  half 
the  computer  time  required  by  the  first  order  gradient 
program.  Intermediate  8(t)  curves  (numbered  by  descent} 
are  shown  in  Fig.  2.  Those  labeled  6  or  less  come  from 
the  penalty  function  stage.  The  F-4  curve  typifies  the 
refinement  stage.  The  converged  curve  is  labeled  F. 

Experiments  were  made  with  fewer  descent  cycles  in 
the  penalty  function  stage.  These  cases  all  required 
more  computer  time  indicating  that  it  is  best  to  be  close 
to  the  desired  extremal  before  going  over  to  the  refine¬ 
ment  stage. 

Generalized  Newton-Raphson  Method1 


The  last  of  the  three  methods  discussed  in  Part  I 
is  based  on  an  algorithm  for  solving  the  two-point  bound¬ 
ary  value  problem  associated  with  the  Euler -Lagrange 
equations.  As  outlined  in  (3),  it  required  a  fixed 
final  time.  Since  the  sample  problem  of  this  paper  is  a 
minimum  time  problem,  the  procedure  was  altered  to  suit 
this  case.  What  follows  is  a  brief  description  of  the 
modified  procedure  and  a  discussion  of  the  numerical  re¬ 
sults  obtained  by  application  of  the  method  to  the  sample 
orbit  transfer  problem. 

The  two-point  boundary  value  problem  resulting 
from  the  Euler-Lagrange  equations  is  given  by 


This  algorithm  was  apparently  first  suggested  for 
boundary  value  problems  by  Hestenes  (j>)  ,  and  was  de¬ 
veloped  further  by  Kalaba  (6_)  .  A  convergence  proof  for 
N  dimensional  systems  was  given  by  McGill  and  Kenneth 
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Control  Angle  Programs  From  Second  Variation  Technique 
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r  =  u 


-  £<« 


a  =  7"  '  JT  + 

r 


A 


u 


! 

2  2-N2 

'  +  v) 


=  f(2) 


C'u 


v  =  - 


~  +  a  (t) 
r  '  ' 


2v  2 


OS + 


-  f<3> 


A  = 
r 


^  A 


2  "  2  3)  Au  2  v 

r  r 


-  f(4) 


•  v 

A  =  -  A  +  -  A 

u  r  r  v 


v  u  , 

A  =  -  2  -  A  +  -  A 
v  r  u  r  v 


-  £<» 
„  f<6> 


where 


a  (t)  =  - t— t— 

iHq  +  mt 


and  the  boundary  conditions  are 
t  =0 


t  =  tj-  (unspecified) 


r(0)  =  r 


0 


r(tf)  =  rf 


u(0)  =  u 


0 


u(tf)  =  uf 


v  (0)  =  V, 


0 


v(tf)  =  vf 


Equations  (9)  through  (14)  may  be  written  as 


X  =  F(X,  t) 


(9) 

(10) 

(11) 

(12) 

(13) 

(14) 

(15) 


(16) 


(17) 
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where 


X  =  (*<«  ,  .  .  . , 

x«) 

F  =  (f(1),  .... 

f  (6)}  <18> 

and 

*(1>  (t)  =  r(t) 

,  x(2)(t)=u(t) 

,  x(3)  (t)  =  v(t) 

x(4)(t)  =  Ar(t) 

,  x(5)(t)  =  Au(c) 

,  X^6'(t)  =  \(.t) 

The  method  proceeds  by  solving  the  following  sequence  of 
linear  two-point  problems. 


n+1 


-  J(v 


t) 


[  Xn+1 


-  X 


nj 


+ 


F<V 


t)  n  =  0,  1, 


(19) 


where  J(X,t)  is  the  Jacobian  matrix  of  partial  deriva¬ 
tives  of  the  f(l)  with  respect  to  the 


x^,  i  =  1,  .....  6,  j  =  1,  .  ..,  6.  A  starting  vector, 
X^(t)  and  an  estimated  final  time,  tf^,  are  assurae<3 

and  the  sequence  of  linear  boundary  value  problems  is 
solved  numerically  by  the  procedure  outlined  in  (3) . 
Altering  the  original  boundary  conditions  so  that  the 

final  r  is  maximized  at  tr  (k=0)  we  have: 

tk 


n 


x 

n 


x 

n 


(4) 
x  v  ' 


t  =  0 

c  =  tfk 

(0)  =  rn(0) 

r0 

(0)  -  un(0) 

o 

d 

ii 

(2) 

X 

n 

(tf)  -  un(tf) 

=  uf 

o 

c 

> 

II 

o 

v0 

(3) 

Xs  ' 

n 

(tf)  =  Vh) 

=  vf 

(0)  =  Ar  (0) 

=  1 

n 


(20) 


n  =  1,  2, 
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Setting  ^(O)  =  1  accomplished  the  scaling  of  the  mul¬ 
tipliers.  The  iteration  proceeds  until  "p(X  , ,  ,  X  )  <  £ 

n-t-i  n  — 

where 


P(X 


n+1 ' 


X  ) 
n' 


I 


max 
te[0,tf  ] 
k 


x 


U)  .  „(i) 


n+1 


n 


(21) 


At  this  stage  the  final  time,  tf  ,  is  adjusted  auto- 

k 

matically  according  to  the  difference  [r£  -  r(t£  )  J  by 

k 

a  scalar  application  of  the  Newton-Raphson  procedure: 


t 


(tfk 

t£k  +  r^f  > 


r(t  ) 
rk-l 


r  (tr  )] 
k 


(22) 


The  above  iteration  on  X^ 

final  time  tf  until  p 
£k+l 


now  continues  for  the  new 
is  again  <  e.  The  over -a 11 


process  proceeds  until  p  <  e  where 


P  P  +  j~|tf 


-  t, 


k+1 


(23) 


and  b  is  a  scaling  factor.  The  corresponding  iterate 
X^^  is  accepted  as  the  solution  to  the  minimum  time 

problem,  and  a  final  check  is  run  by  integrating  the  non¬ 
linear  Euler- Lagrange  equations  with  a  complete  set  of 
initial  conditions  taken  from  the  final  iterate. 

For  the  purpose  of  numerical  precision  the  data 
for  the  sample  problem  were  normalized  to  obtain 


1.000 

v£  = 

.8098 

1.523 

uf  = 

0.000 

1.000 

in0  = 

1.000 

1.000 

m  = 

-  .07487 

0.000 

T  = 

.1405 

(24) 
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This  resulted  in  a  time  unit  of  58.18  days.  The  starting 
vector  X{)(t)  was  chosen  rather  crudely  as  follows: 


t 


-  178.0  days,  or  3.0G0  of  our  time  units 


xQ(1)  (t)  =  rQ(t) 


rf  '  ro 
°  ro  +  T 


Xq2^  (t)  =  uQ  (t)  =  0 

x0(3)(t)  =  v0(t)  = 

xi4)(t)  =  A  (t)  =  1.000 

0  '  r0 


(25) 


<o 


(t) 


j  .5200 
l-. 5000 

j  .3000 

VoOO 


for  te [0 ,  i  tf  J 
for  te[|  tfQ,  tfQ] 

for  te[0,  j  tf  J 

for  tfe'ti  tf  ,  tf  ] 
0  0 


The  final  two  starting  functions  \,  (t)  and  K.  (t) 

0  0 

correspond  to  a  control  angle  vq  (t)  which  is  constant 
at  60  above  the  local  horizontal  for  the  first  half  of 
the  transit  time,  and  constant  inward  along  the  local 
vertical  for  the  remaining  half  of  the  transit  time  (see 
Fig.  3). 

The  sequence  (X^  converged  uniformly  to  an  ac¬ 
curacy  of  5  significant  figures  in  13  total  iterations 
with  4  shifts  of  the  final  time.  The  resultant  minimum 
time  was  found  to  be  193.2  days;  in  agreement  with  the 
results  of  the  gradient  and  second  variation  solutions. 
The  total  computer  time  (IBM  '  094)  required  was  36  sec¬ 
onds.  Figure  3  illustrates  the  behavior  of  the  control 
angle  program,  where  eQ(t)  *-£  the  starting  function. 
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l (t)  through  4(t)  correspond  to  the  4  shifts  of  the 
final  time  tf,  and  ,v(t)  results  from  the  integration 

of  the  nonlinear  Euler-Lagrange  equations  with  the  ini¬ 
tial  values  taken  from  the  final  iterate.  The  curves 
for  2(0*  3(f)  ,  and  ^(t)  lie,  within  our  plotting 

accuracy,  on  the  solution  curve  t  “(t)  ;  except  for  the 
final  segments  as  indicated  on  the  figure.  Figure  4 
illustrates  the  behavior  of  the  metrics,  p,  and  ,  . 

We  observe  that  for  this  particular  example  the 
approach  just  described  is  systematic,  simple  to  apply, 
and  yields  rapid  convergence  from  crude  a  priori  starting 
functions.^-  However,  it  is  possible  that  for  other  var¬ 
iational  problems  a  priori  starting  functions  sufficient 
to  produce  convergence  may  not  be  easily  obtainable.  In 
such  a  case  one  might  consider  the  possibility  of  using 
a  hybrid  approach,  that  is,  a  combination  of  one  of  the 
gradient  procedures,  or  a  crude  application  of  dynamic 
programming;  and  the  technique  just  described.  Such  a 
hybrid  approach  is  currently  under  study. 

Acknowledgments 


The  work  reported  in  this  paper  was  partly  sup¬ 
ported  by  the  USAF  Office  of  Scientific  Research,  Applied 
Mathematics  Division,  under  Contract  AF  49 (638) -1207  and 
by  NASA  Marshall  Space  Flight  Center,  Aero-Astrodynamics 
Laboratory,  under  Contract  NAS  8-1549. 

The  authors  wish  to  express  their  appreciation  to 
Dr.  Henry  J.  Kelley  for  valuable  discussions  and  sugges¬ 
tions  concerning  this  work;  and  to  Gerald  E.  Taylor 
whose  able  programming  made  the  numerical  example  for 
the  generalized  Newton-Raphson  method  possible. 

References 


of 


(1)  .  Lindorfer,  W.,  and  Moyer,  H.G.,  "Application 
Low  Thrust  Trajectory  Optimization  Scheme  to 
Planar  Earth-Mars  Transfer,"  ARS  Journal ,  February 
1962. 


For  additional  numerical  examples  see  (_Z)  . 

10 


<•.***■ 


COMPUTING  METHODS  IN  OPTIMIZATION  PROBLEMS 


(2J).  Kelley,  H.J.,  Kopp,  R.E.,  and  Moyer,  H.G.,  "A  Tra¬ 
jectory  Optimization  Technique  Based  upon  the 
Theory  of  the  Second  Variation, "  presented  at  the 
AIAA  Astrodynamics  Conference,  Yale  University, 
Connecticut,  August  1963. 

(3) .  McGill,  R.,  and  Kenneth,  P.,  "A  Convergence  Theorem 

on  the  Iterative  Solution  of  Nonlinear  Two-Point 
Boundary  Value  Systems,"  presented  at  the  XIVth 
IAF  Congress,  Paris,  France,  September  1963. 

(4) .  Kelley,  H.J.,  "Method  of  Gradients,"  Chapter  6  of 

Optimization  Techniques,  edited  by  G.  Leitmann, 
Academic  Press,  1962. 

(5) .  Hestenes,  M.R.,  Numerical  Methods  of  Obtaining  So¬ 

lutions  of  Fixed  End  Point  Problems  in  the  Calculus 
of  Variations.  RM-102,  The  Rand  Corporation,  August 
1949.  . 

(6)  .  Kalaba,  R.,  "On  Nonlinear  Differential  Equations, 

the  Maximum  Operation,  and  Monotone  Convergence," 
Journal  of  Mathematics  and  Mechanics,  Vol.  8, 

No.  4,  pp.  519-574,  July  1959. 

(JJ.  McGill,  R.,  and  Kenneth,  P.,  "Solution  of  Varia¬ 
tional  Problems  by  Means  of  a  Generalized  Newton- 
Raphson  Operator,"  to  be  published. 


11 


